function M=mut_ind(h)
%h is usually negative
Nc=1;
Np=10;
rc=1;
rp=0.9;
lc=0.1;
lp=10;

%calculate the double integral:
function out = integrnd(t1,t2)
    out =  rc*rp*cos(t1-t2)./sqrt(rc^2+rp^2-2*rc*rp*cos(t1-t2)+(h+lc*t1/2/pi/Nc+lp*t2/2/pi/Np).^2);
end

M = 1e-7*dblquad(@integrnd,-pi*Nc,pi*Nc,-pi*Np,pi*Np,1e-2);

end